# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <errno.h>
# include <math.h>
/*
# ifdef _OMP_SOLVE_
# include <omp.h>
# endif*/

# ifdef SOLVE_ENERGY
extern double temp_advance(struct grid *g, const struct arr_struct *depVar);
# endif

extern void init_cond(struct grid *g);
extern short print_result(unsigned int t_cnt, const struct grid *g);
extern unsigned int vel_advance(struct grid *g, unsigned int *pcnt_max, const struct arr_struct *depVar, const struct arr_struct *memFace);
extern double inlet_velocity(double r, double t, short index);
extern double dens(double T);
extern double kvisc(double T);
extern double thm_cond(double T);
extern short print_recov(const struct grid *g);
extern double orlansky_outlet(struct grid *g, unsigned int t_cnt, const struct arr_struct *memOut);

static double heat_balance(const struct grid *g, double E_absr);
static FILE *file_open(FILE *lptr, const char *file, unsigned int t_cnt, unsigned int nfreq);


extern short solve(struct grid *g, const struct arr_struct *depVar, const struct arr_struct *memFace, const struct arr_struct *memOut)
 {
  unsigned int t_cnt, pcnt_max, pcnt_total;
  short sst_cnt, div_cnt;
  double pm_accum, pE_accum, E_absr;

# ifdef TRACE_FILE
  fprintf(g->trptr, "Entering into 'solve()'...\n");
# endif

# ifdef SOLVE_ENERGY
  g->vhsrc_flag = VHSRC_FLAG;
# endif

  g->density0 = dens(T_INLET);
  g->kvisc0 = kvisc(T_INLET);
  g->Re0 = 0.5*W_MAX*d_NOZZLE/g->kvisc0;
  fprintf(g->lptr, "\nReynolds Number at Inlet, Re0 = %E\n", g->Re0); fflush(g->lptr);
  g->thmcond0 = thm_cond(T_INLET);

  init_cond(g);

/* =================== Calculation at different time-steps ==================== */

  pcnt_total = sst_cnt = div_cnt = 0;
  g->dptr = g->tptr = (FILE *)0;
  t_cnt = 0;

  while(g->t <= tMAX)
    {
     if(t_cnt%DFILE_FREQ == 0)
          {
           if(g->dptr) if(g->dptr != stderr) fclose(g->dptr);
           g->dptr = file_open(g->lptr, DISPLAY_FILE, t_cnt, DFILE_FREQ);
          }
     fprintf(g->dptr, "\n\n\t[[ t_cnt = %u, t = %E Sec ]]\n", t_cnt, g->t); fflush(g->dptr);

     if(t_cnt%IRFILE_FREQ == 0)
          {
           if(g->tptr) if(g->tptr != stderr) fclose(g->tptr);
           g->tptr = file_open(g->lptr, ittr_RES_FILE, t_cnt, IRFILE_FREQ);
          }
     fprintf(g->tptr, "\n\n[[ t_cnt = %u, t = %E Sec ]]; ", t_cnt, g->t); fflush(g->tptr);

     pm_accum = orlansky_outlet(g, t_cnt, memOut);
     fprintf(g->tptr, "pm_accum = %E%c; ", pm_accum, '%'); fflush(g->tptr);

     pcnt_total += vel_advance(g, &pcnt_max, depVar, memFace);

#    ifdef SOLVE_ENERGY
     E_absr = temp_advance(g, depVar);
     pE_accum = heat_balance((const struct grid *)g, E_absr);
     fprintf(g->tptr, "pE_accum = %E%c; ", pE_accum, '%'); fflush(g->tptr);
#    endif

     if(pcnt_max == POISS_CNT_MIN) ++sst_cnt; else sst_cnt = 0;
     if(pcnt_max == POISS_CNT_MAX+1) ++div_cnt; else div_cnt = 0;

     t_cnt++;
     g->t += g->dt;

  /* Printing result and recovery files. */
     if( t_cnt%RESF_FREQ == 0 || sst_cnt == SST_CNT_MAX ) print_result(t_cnt, g);

     if(pcnt_total >= RECOVF_POISS_CNT) if(pcnt_max <= POISS_CNT_MAX)
          {
           print_recov(g);
           pcnt_total = 0;
          }

     if(sst_cnt == SST_CNT_MAX)
           {
            fprintf(g->tptr, "STEADY-STATE solution; "); fflush(g->tptr);
            fprintf(g->lptr, "\n\t[Solution reached steady-state]\n"); fflush(g->lptr);
            break;
           }

     if(div_cnt == DIV_CNT_MAX)
           {
            fprintf(g->tptr, "DIVERGED solution; "); fflush(g->tptr);
            fprintf(g->lptr, "\n\t[Solution diverged\n"); fflush(g->lptr);
            break;
           }
    }
/* --------------------------------------Time-loop ends------------------------------------------ */

  if(g->dptr) if(g->dptr != stderr) fclose(g->dptr);
  if(g->tptr) if(g->tptr != stderr) fclose(g->tptr);

# ifdef TRACE_FILE
  fprintf(g->trptr, "...\"solve()\" ends.\n\n");
# endif

  return g->t > tMAX ? 1:-1;
 }



# ifdef SOLVE_ENERGY
static double heat_balance(const struct grid *g, double E_absr)
    {
     double E_inlet, E_outlet, E_vhsrc, E_heater, *th, *r;
     unsigned short i, j, ii, jj, nth;

#    ifdef TRACE_FILE
     fprintf(g->trptr, "Entering into 'heat_balance()'...\n");
#    endif

     nth = g->nth; th = g->th; r = g->r;

  /* Determining total energy created by the volume-heat-source. */
     E_vhsrc = (g->vhsrc_flag == 'y' ? \
                  (th[NthG+nth]-th[NthG])*d_CYL*d_CYL*h_CYL*g->density0*Cp_WATER*(EXPECTED_T_RISE/EXPECTED_TIME)*g->dt/8:0);

  /* Calculating heat input because of heater on the top surface. */
     E_heater = g->dt*TOP_HEATER_HF*0.5*g->r_heater*g->r_heater*(g->th[NthG+g->nth]-g->th[NthG]);

  /* Calculating the energy coming into the domain through the inlet in time t~t+dt. */
     E_inlet = 0.25*Cp_WATER*T_INLET*g->density0*r_NOZZLE*r_NOZZLE*W_MAX\
       *(g->dt-tFRAC*(exp(log(1-wFRAC)*(g->t+g->dt)/tFRAC) - exp(log(1-wFRAC)*g->t/tFRAC))/log(1-wFRAC))*(th[NthG+g->nth]-th[NthG]);

  /* Calculating the heat gone out of the domain through outlet between t~t+dt. */
     for(E_outlet = 0, i = 0, ii = NrG + g->nr - g->nr_gap; i < g->nr_gap; i++, ii++) for(j = 0, jj = NthG; j < g->nth; j++, jj++)
     E_outlet += -0.0625*(r[ii+1]+r[ii])*(r[ii+1]-r[ii])*(th[jj+1]-th[jj])*(g->outBC.w[0][i][j][0]+g->outBC.tnext_wg[i][j])\
                     *(g->outBC.T[0][i][j][0]+g->outBC.T[0][i][j][1]+g->outBC.tnext_Tg[i][j]+g->T[ii][jj][NzG]);
     E_outlet *= g->density0*Cp_WATER*g->dt;

#    ifdef TRACE_FILE
     fprintf(g->trptr, "...'heat_balance()' ends.\n");
#    endif

  /* Return % heat-accumulation in the domain between time t~t+dt. */
     return (1-(E_outlet+E_absr)/(E_inlet+E_heater+E_vhsrc))*100;
    }
# endif


static FILE *file_open(FILE *lptr, const char *file, unsigned int t_cnt, unsigned int nfreq)
    {
     FILE *ptr;
     char fname[MAX_FNAME_LEN], err_str[SYSERR_STR_LEN];

     if(file == (char *)0) ptr = stderr;
     else{
          if(!strcmp(file, NULL_DEV)) strcpy(fname, NULL_DEV);
          else sprintf(fname, "%02d.%s.%09d-%09d.dat", RUN_CNT, file, t_cnt, t_cnt+nfreq-1);

          if( (ptr = (FILE *)fopen(fname, "w")) == (FILE *)0)
                {
                 strerror_r(errno, err_str, SYSERR_STR_LEN);
                 fprintf(lptr, "WARNING:: In 'file_open()': Unable to open '%s' for writing: %s.\nRedirecting output to 'stderr'.\n",\
                                        fname, err_str);
                 ptr = stderr;
                }
         }

     return ptr;
    }
